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1  Task  Objects 

Drag  reduction  using  polymer  additives  in  wall-bounded  flows  poses  many  challenges  to  our 
understanding  of  turbulence  and  polymer  dynamics.  Up  to  80%  of  friction-drag  reduction 
can  be  obtained  with  ultra  dilute  solutions  of  high  molecular  weight  polymers,  for  which 
each  molecule  is  several  orders  of  magnitude  smaller  than  the  smallest  turbulent  scale  of 
the  flow.  To  investigate  this  multiscale  problem,  numerical  simulations  at  the  microscopic 
and  macroscopic  level  were  performed.  At  the  microscopic  level,  Brownian  dynamics  (BD) 
simulations  were  used  to  assess  the  realism  of  polymer  models  and  to  document  the  effect 
of  turbulence  on  polymer  configuration.  Drag  reduced  flows  were  investigated  using  direct 
numerical  simulations  (DNS)  of  viscoelastic  flow  (macroscopic  simulation).  DNS  was  able 
to  show  how  energy  is  transferred  between  turbulent  structures  and  polymers  in  turbulent 
channel  flow  and  boundary  layers.  It  also  allowed  for  the  detailed  investigation  of  the  two 
different  regimes  of  polymer  drag  reduction,  namely  low  drag  reduction  (LDR)  and  high 
drag  reduction  (HDR).  Based  on  the  DNS  and  BD  results,  a  model  for  the  mechanism  of 
polymer  drag  reduction  has  been  proposed.  The  I)NS  formalism  was  also  extended  to  fibres 
additives  for  which  a  model  of  the  mechanism  of  drag  reduction  was  also  derived. 


2  Technical  problem 

BD  and  DNS  were  solved  using  state  of  the  art  flow  solvers  developed  at  Stanford.  The 
challenge  of  the  BD  simulation  was  to  implement  an  efficient  and  reliable  particle  tracking 
technique  in  a  fully  parallelized  code.  The  DNS  of  viscoelastic  flow  had  to  overcome  small- 
scale  instabilities  arising  from  the  resolution  of  the  macroscopic  polymer  model.  Both  of 
these  problems  were  solved  adequately  as  is  described  in  the  Technical  Results  section. 

3  General  Methodology 

The  addition  of  small  amounts  of  long  chain  polymer  molecules  to  wall-bounded  flows  can 
lead  to  dramatic  drag  reduction.  Although  this  phenomenon  has  been  known  for  about  fifty 
years,  the  action  of  the  polymers  and  its  effect  on  turbulent  structures  are  still  unclear. 
Detailed  experiments  have  characterized  two  distinct  regimes  (Warholic  et  al,  1999),  which 
are  referred  to  as  low  drag  reduction  (LDR)  and  high  drag  reduction  (HDR).  The  first  regime 
exhibits  similar  statistical  trends  as  Newtonian  flow:  the  log-law  region  of  the  mean  velocity 
profile  remains  parallel  to  that  of  the  Newtonian  flow  but  its  lower  bound  moves  away  from 
the  wall  and  the  upward  shift  of  the  log-region  is  a  function  of  drag  reduction,  herein  referred 
to  as  DR.  Although  streamwise  fluctuations  are  increased  and  transverse  ones  are  reduced, 
the  shape  of  the  rms  velocity  profiles  is  similar.  At  higher  drag  reductions,  larger  than  about 
40%,  the  flow  enters  the  HDR  regime  for  which  the  slope  of  the  log-law  is  dramatically 
augmented  and  the  Reynolds  shear  stress  is  small  (Warholic  et  a/.,  1999;  Ptasinski  et  aZ., 
2003).  The  drag  reduction  is  eventually  bounded  by  a  maximum  drag  reduction  (MDR,  Virk 
&  Mickley,  1970)  which  is  a  function  of  the  Reynolds  number. 

In  drag  reduced  flows,  a  stress  deficit  is  observed  in  the  stress  balance  whose  large  magni¬ 
tude  at  HDR  has  been  interpreted  as  the  necessary  input  of  energy  from  the  polymers  to  the 
flow  for  the  sustenance  of  the  asymptotic  MDR  turbulence  (Warholic  et  aZ.,  1999).  Recently, 
numerical  simulations  have  allowed  the  simultaneous  study  of  velocity  and  polymer  fields, 
giving  the  opportunity  to  relate  turbulence  and  polymer  stress  directly.  Several  explanations 
have  been  proposed  for  the  mechanism  of  polymer  drag  reduction  through  such  computer 
tions.  Based  on  averaged  budgets  of  kinetic  and  energy,  Dimitropoulos  et  al  (2001)  observed 
that  streamwise  enstrophy  is  inhibited  by  the  extensional  viscosity  generated  by  polymers 
stretched  by  turbulence.  Based  on  a  similar  argument,  Ptasinski  et  al  (2003)  related  the 
damping  of  vortices  to  a  shear  sheltering  effect  occurring  in  the  near-wall  region  which  would 
effectively  decouple  the  outer  layer  vortices  from  the  inner  layer.  Neither  paper  proposes  a 
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mechanism  for  the  injection  of  energy  from  polymers  to  turbulence.  Min  et  al.  (20036, a) 
analyzed  their  simulation  in  terms  of  elastic  energy,  observing  that  a  significant  transport  of 
this  energy  from  the  viscous' sublayer  to  the  buffer  layer  and  the  log  region.  They  derived  a 
criterion  of  the  onset  of  drag  reduction  based  on  the  relaxation  time  and  the  time  required 
for  polymers  to  travel  from  the  viscous  sublayer  to  the  buffer  layer.  These  authors  explain 
MDR  by  the  release  of  energy  into  the  buffer  layer,  which  would  be  sufficient  to  sustain 
turbulence. 

All  the  studies  referenced  above  have  relied  on  DNS  of  turbulent  channel  flow  and  conclu¬ 
sions  were  derived  from  time  and  space  averaged  statistics.  Such  statistics  can  hide  important 
intermittent  events.  Due  to  the  intermittent  character  of  near-wall  turbulence,  our  effort  fo¬ 
cused  on  the  accurate  computation  of  small  scale  events  and  their  characterizations  in  terms 
of  coherence  and  relation  to  turbulent  structures.  BD  using  advanced  polymer  models  was 
used  for  the  first  time  in  turbulent  channel  flows,  providing  insights  on  the  dynamics,  while 
DNS  was  able  to  document  local  transfer  of  energy  between  polymers  and  turbulence.  This 
combined  effort  led  to  a  model  of  polymer  drag  reduction  which  describes  how  turbulent 
energy  is  stored  by  polymers  and  where  this  energy  is  released,  for  both  LDR  and  HDR 
regimes.  Simulations  of  turbulent  boundary  layers  were  developed  for  comparison  with  the 
Stanford  experiment  in  Phase  II  of  the  FViction  Drag  Reduction  project.  Finally,  DNS  of 
fibre  additives  highlighted  a  different  mechanism  of  drag  reduction. 

3.1  Polymer  models 

All  the  models  tested  here  are  derived  from  the  Finitely  EJxtensible  Non-linear  Elastic 
(FENE)  model.  Different  levels  of  refinement  were  tested  using  BD  simulations,  while  DNS 
was  performed  using  the  FENE-P  model  only. 

3.1.1  FENE 

The  FENE  model  is^  based  on  an  elastic  dumbbell  which  represents  the  polymer  molecule 
as  two  beads  connected  by  a  spring  Bird  et  al.  (1987);  Herrchen  &  Oettinger  (1997).  The 
non-linear  spring  force  is  modeled  by  a  Warner  spring  representing  the  entropic  forces  which 
tend  to  bring  the  polymer  back  into  its  coil  configuration.  The  non-linearity  of  the  spring 
ensures  that  the  dumbbell  cannot  extend  more  than  a  given  maximum  extensibility.  The 
equation  for  the  extension  of  the  polymer  is  derived  by  neglecting  the  inertia  of  the  beads  and 
equilibrating  the  different  forces  acting  on  the  beads.  Polymer  molecules  (microscopic  model) 
or  solutions  (macroscopic  or  constitutive  model)  axe  defined  by  the  polymer  relaxation  time 
(Weissenberg  number  We)  and  the  maximum  length  of  a  polymer  molecule  L.  In  the  case 
of  a  polymer  solution,  /3,  the  ratio  of  the  solvent  viscosity  v  to  the  total  viscosity,  effectively 
controls  the  concentration  of  polymers. 

3.1.2  FENE-P 

The  FENE-P  model  (P  for  Perterlin)  describes  the  ensemble  average  evolution  of  polymer 
configuration,  assuming  homogeneous  concentration.  The  FENE-P  model  is  therefore  macro¬ 
scopic  and  generally  described  in  the  Eulerian  framework.  BD  simulations  using  this  model 
were  also  performed  to  compare  Eulerian  and  Lagrangian  advection  of  polymer  quantities. 

3.1.3  FENE-P  Brownian 

This  model  is  hybrid  FENE/FENE-P  model,  where  the  average  polymer  configuration  is 
computed  using  FENE-P  but  the  model  still  has  the  stochastic  Brownian  term. 

3.1.4  FENE  MULTICHAIN 

The  FENE  multichain  model  It  consists  of  a  chain  of  beads  connected  by  springs.  The 
derivation  of  the  equation  for  the  end-to-end  vectors  is  similar  to  that  of  the  FENE  dumbbell 
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Somasi  et  al  (2002a).  The  bead-spring  chain  is  a  less  coarse  approximation  of  a  polymer 
molecule,  since  it  takes  account  of  some  of  the  internal  relaxation  modes  which  are  present 
in  a  real  molecule. 

3.1,5  Fibres  model 

Drag  reduction  induced  by  rigid  fibres  in  a  turbulent  channel  flow  using  direct  numerical  sim¬ 
ulation  were  investigated  using  the  channel  fiow  code  developed  for  polymer  drag  reduction. 
The  fibres  are  modeled  using  a  well-known  constitutive  equation  involving  the  moments  of 
orientation  vector  and  two  closure  approximations  (Hinch  &  Leal,  1976).  The  model  has  a 
form  similar  to  the  constitutive  FENE-P  model. 

3.2  Flow  solver 

The  fiow  solver  is  based  on  robust  finite  difference  schemes  which  allow  for  the  simulation  of 
complex  geometries.  The  code,  based  on  Pierce  (2001),  utilizes  a  staggered  grid  and  second- 
order  finite  differences  for  spatial  derivatives.  Time  integration  is  semi-implicit,  using  either 
a  third  order  Runge-Kutta/  Crank  Nicolson  scheme  Dubief  (2002)  or  an  iterative  method 
developed  by  Pierce  (2001)  to  enhance  secondary  conservation  and  nonlinear  stability.  Sim¬ 
ilar  codes  have  been  extensively  used  and  validated  at  the  Center  for  Turbulence  Research 
for  a  wide  variety  of  problems  including  combustion  and  particle-laden  flows. 

3.3  Particle  tracking 

For  BD  simulations,  polymer  trajectories  are  solved  using  Lagrangian  advection.  The  in¬ 
tegration  of  the  polymer  molecule  trajectories  assumes  that  the  molecule  center-of-mass 
motion  is  characterized  by  no  inertia  and  an  infinite  Peclet  number.  Therefore,  these  tra¬ 
jectories  represent  an  exact  Lagrangian  description  of  the  flow.  To  accommodate  the  small 
scale  behavior  of  polymers,  the  trajectories  are  integrated  at  a  smaller  time  step  than  that 
of  the  flow.  Therefore,  the  velocity  at  the  intermediate  time  steps  is  obtained  by  linearly 
interpolating  two  flow  fields  in  time.  The  polymer  molecule  itself  is  advanced  using  a  second 
order  Runge-Kutta  scheme.  Since  the  Eulerian  velocity  is  only  known  on  a  grid  network, 
a  spatial  interpolation  is  needed,  to  determine  the  velocity  at  its  current  position.  The  tri- 
linear  interpolation  was  found  to  be  sufficient  (Terrapon  et  aZ.,  2003).  The  time  integration 
of  FENE  models  uses  a  semi-implicit  predictor-corrector  scheme  (Oettinger,  1996;  Somasi 
et  aLj  20026) 

3.4  DNS  of  viscoelastic  flows 

The  FENE-P  model  in  the  Eulerian  framework  leads  to  the  calculation  of  polymer  stress 
which  enters  the  Navier  Stokes  equations  as  an  extra  term.  The  numerical  solution  of  the 
FENE-P  model  requires  high  order  methods  as  pointed  out  by  Min  et  al,  (2001).  To  capture 
the  small  scales  generated  by  the  model,  some  derivatives  have  to  be  calculated  using  compact 
schemes  Dubief  (2002),  as  is  discussed  in  the  Technical  Results  section. 

3.5  Flow  description 

3.5.1  Channel  flow 

The  channel  flow  domain  assumes  periodicity  in  directions  parallel  to  walls,  streamwise  and 
spanwise.  No-slip  is  prescribed  on  walls.  Mass  flow  is  kept  constant  so  that  drag  reduction 
is  measured  as  a  drop  in  streamwise  pressure  gradient. 

DNS  simulations  are  performed  at  an  intermediate  Reynolds  number,  Rcc  =  7500,  or  a 
bulk  Reynolds  number  of  ReM  =  5000.  Conservation  of  the  mass  flow  is  imposed  which  gives 
=  hUrjv  =  300,  where  h  is  the  channel  half  width,  Ur  =  \/^{dUxldy)^aii  is  the  friction 
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velocity  and  u  is  the  kinematic  viscosity  The  superscript  indicates  quantities  normalized 
by  Ur  and  u.  The  resolution  is  Aa?"*"  =  9,  ==  0.1  —  5  and  =  6,  when  normalized  by 

the  skin  friction  at  DR=0%.  The  polymer  parameters  used  in  the  course  of  this  study  are 
shown  in  Table  1.  The  streamwise  and  spanwise  dimensions  and  resolutions  were  adjusted 
to  accommodate  highly  drag  reduced  flows  which  length  scales  grow  significantly.  According 
to  the  Stanford  experiment,  the  streamwise  coherence  of  the  streaks  could  be  increased  by 
an  order  of  magnitude  of  the  Newtonian  length  scale  1000.  In  a  preliminary  parametric 
study  for  60%  of  drag  reduction,  streamwise  correlations  of  velocity  were  indeed  found  to  be 
small  yet  finite  even  for  Lj  =  6000  based  on  the  Newtonian  friction  velocity,  while  =  1000 
proved  to  be  sufficient.  Decreasing  to  4000  had  a  negligible  impact  on  the  magnitude  of 
drag  reduction,  less  that  5%  well  within  the  convergence  error  of  statistics,  and  this  length 
was  adopted  for  the  present  study.  The  dimensions  of  the  computational  domain  is  thus 
47rft  X  2h  X  4h, 

For  BD  simulations,  the  computational  domain  is.  reduced  to  a  minimal  channel  flow 
(Jimenez  &  Moin,  1991)  of  dimensions  7thx2hx  h.  The  simulation  has  been  carried  out  on 
a  grid  of  96  x  151  x  64  points  so  that  dx'^  =  9.5,  =  0.2  -  10  and  =  4.5.  Finally  a 

constant  time  step  of  dt  =  0,01H/U  was  chosen,  ensuring  a  maximum  CFL  =  0.6.  To  achieve 
converged  statistics,  a  large  ensemble  of  particles  are  needed.  In  this  case  N  =  10®  particles 
were  used  with  each  of  them  having  a  different  trajectory.  The  extensibility  parameter  was 
chosen  to  be  =  L^b  =  3600  in  order  to  simulate  real  molecules.  Calculations  were  also 
performed  for  b  =  900  and  b  =  10000,  which  gave  qualitatively  the  same  results.  The  bead¬ 
spring  chain  used  had  5  springs  corresponding  to  5  relaxation  modes.  In  all  simulations  we 
used  a  Weissenberg  number  We  =  3  corresponding  to  a  Weissenberg  number  based  on  the 
wall  shear  We~^  =  =  34.  The  equations  for  the  polymer  molecules  were  advanced 

using  a  time  step  dtp  =  0.001  which  was  ten  times  smaller  than  that  of  the  flow  calculation. 

3.5.2  Boundary  layer 

For  the  momentum  equations,  the  inflow,  outflow  and  fi:ee-stream  boundary  conditions  are 
very  important  given  that  they  are  approximate  and  not  trivial.  For  the  outflow,  the  convec¬ 
tive  condition  is  used  (see  the  discussion  in  Na,  1996).  For  the  free-stream,  the  zero  vorticity 
condition  is  applied  for  the  streamwise  velocity,  w,  the  zero  shear  constraint  for  the  spanwise 
velocity,  w  and  the  outflow  velocity,  v,  is  consistent  with  the  boundary  layer  growth,  through 
on  the  fly  computation  of  an  integral  scale  (displacement  thickness)  (Na,  1996;  Lund  et  a/., 
1998;  Jacobs  &  Durbin,  2001): 

No  slip  and  periodic  conditions  are  used  in  the  spanwise  direction  and  on  the  wall. 
The  inflow  condition  is  supplied  from  a  Newtonian  simulation  which  uses  an  inflow  profile 
constructed  by  recycling  from  a  downstream  station,  using  the  approach  of  Lund  et  al. 
(1998).  The  recycling  method  of  Lund  has  been  proven  to  minimize  errors  at  the  inflow 
due  to  flow  readjustment.  Alternatively,  one  can  use  recycling  in  conjunction  with  a  freely 
developing  inflow  conformation  profile. 

Several  simulations  of  viscoelastic  turbulent  boundary  layer  flow  with  the  FENE)-P  model 
have  been  performed  recently  in  order  to  validate  the  method.  The  first  simulations  were 
for  uncoupled  flow  at  a  maximum  extension  of  900  and  a  Weissenberg  number  (based  on 
the  Newtonian  inflow  shear  rate  at  the  wall),  We^  equal  to  25.  The  coupled  simulations 
were  first  at  the  extensibility  and  elasticity  values  described  above  and  then  at  =  3600 
and  Wer  =  50.  The  Reynolds  number,  Re^,  considered  so  far  is  relatively  low,  equal  to  300 
when  based  on  the  inlet  momentum  thickness.  This  was  motivated  by  the  desire  to  minimize 
the  computational  cost  of  the  first  cases  studied  in  order  to  facilitate  the  development  and 
validation  of  the  algorithm.  At  these  parameter  values  a  grid  of  128  x  65  x  128  in  the 
streamwise,  x,  wall-normal,  y,  and  spanwise,  direction  is  adequate  (Na,  1996).  The  size  of 
the  computational  domain  was  equal  to  87^*^  x  19<5*^  x  25 Jin,  where  Jin  is  the  displacement 
thickness  at  the  inflow  boundary.  The  spatial  resolution  in  wall  units  (calculated  with  the 
inlet  friction  velocity)  is  in  all  cases  Ax+  =  18.4,  Ay+  =  (0.17  -  27)  (a  hyperbolic  tangent 
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stretching  function  is  used)  and  Az"*"  =  5.3. 

I 

3.6  Technical  results 

3.7  Brownian  dynamics 

For  this  project,  our  first  approach  was  to  compare  the  different  polymer  models.  Figure  1 
shows  the  mean  and  RMS  of  the  square  extension  non-dimensionalized  by  the  maximum 
extensibility  parameter  6  as  a  function  of  the  distance  to  the  wall  The  square  exten¬ 
sion  of  the  polymer  is  defined  as  the  norm  square  of  the  extension  vector  Q.  For  the 
constitutive  FENE)-P  model,  it  is  defined  as  the  trace  (QiQi)  of  the  conformation  tensor. 
It  can  be  seen  that  all  models  give  qualitatively  the  same  prohle.  However,  the  FENE  and 
MULTICHAIN  models  have  a  slightly  lower  mean  extension.  The  RMS  are  similar  for  all 
models  except  for  the  Brownian  FENE-P,  which  has  much  higher  fluctuations.  A  compari¬ 
son  of  the  different  models  along  the  same  trajectory  is  presented  in  Fig.  2  (note  that  the 
initial  conditions  were  different).  The  behavior  of  the  different  models  is  again  qualitatively 
similar,  even  if  each  particle  has  a  different  stretch  history. 

To  investigate  the  dynamics  of  polymer  in  turbulence,  the  evolution  of  polymer  stretch 
was  found  to  correlate  to  that  of  eigenvalues  and  invariants  of  the  velocity  gradient  (Chong 
et  alj  1990).  For  an  incompressible  flow  the  eigenvalues,  (7,  are  obtained  as  solutions  of  the 
characteristic  equation 

-f  Qa  -h  i?  =  0,  (1) 

with  the  tensor  invariants  Q  and  R  given  by 

Q  =  -l/2tr((Vu)2),  (2) 

R  =  -det(Vu).  '  (3) 

The  nature  of  the  eigenvalues  is  determined  by  the  discriminant  D  =  (27/4)Ji^  D  >  0 
gives  rise  to  one  real,  two  complex-conjugate  eigenvalues;  D  <  0  gives  three  real  distinct 
eigenvalues  and  D  =  0  corresponds  to  three  real  eigenvalues  of  which  two  are  equal. 

Following  this  approach,  it  was  shown  that  regions  of  the  flow  where  D  >  0  and  i?  >  0  in  a 
{Q,R)  plot  (implying  a  negative  real  eigenvalue  and  two  complex-conjugate  eigenvalues  with 
positive  real  parts)  unravel  the  polymer  molecule  at  a  sufficient  high  Weissenberg  number. 
In  this  particular  case,  the  stretching  does  not  occur  along  a  specific  direction  but  along  a 
rotating  axis  in  the  plane  associated  with  the  complex-conjugate  eigenvalues.  This  shows  that 
the  stretching  of  the  polymer  is  driven  by  the  extensional  character  of  the  flow  quantified  by 
the  positive  real  part  of  the  eigenvalues.  From  this  topological  analysis,  one  can  conclude  that 
a  polymer  molecule  will  only  fully  unravel  if  it  experiences  a  flow  with  a  strong  extensional 
character,  i.e.,  one  of  the  eigenvalues  of  the  local  velocity  gradient  tensor  has  a  large  positive 
real  part.  Therefore,  we  introduce 

a*  =  max  (Re(c7i)) .  (4) 

as  a  measure  of  the  ability  of  the  fiow  to  stretch  the  polymer  molecule. 

In  Fig.  3  shows  the  history  of  cr*,  Q  and  R  that  particles  experiences  before  they  reach 
significant  extensions.  The  plots  are  obtained  by  conditional  sampling  based  on  the  extension 
of  polymers.  At  low  We  the  polymer  molecules  experience  a  burst  in  a*  just  before  reaching 
the  threshold  value.  When  We  is  increased,  the  maximum  value  of  a*  decreases  (see  figure 
3  a).  Figure  3(d)  illustrates  that  even  at  a  larger  We  a  large  cr*  is  needed  to  extend  the 
polymer  near  to  its  maximum  extension. 

Flow  visualizations  also  demonstrated  that  the  regions  of  large  cr*  are  always  located 
next  to  the  vortices  and  can  also  be  seen  as  structures  advected  by  the  mean  flow.  This  is 
illustrated  in  figure  4,  which  shows  isosurfaces  of  Q  representing  the  vortices,  isosurfaces  of 
a*  and  the  polymer  molecules  which  are  highly  stretched.  The  correlation  between  a*  and 
the  vortices  is  striking  and  can  provide  new  insight  into  the  mechanisms  of  polymer  drag 
reduction. 
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3.8  DNS  of  channel  flow 

The  stability  of  the  constitutive  equations  was  carefully  investigated  in  Dubief  et  al  (2004a). 
Due  to  the  lack  of  diffusion  in  the  constitutive  FENE-P  model,  the  polymer  configuration 
field  tends  to  develop  sharp  gradient.  As  a  consequence,  an  approach  similar  to  that  used 
for  shocks  is  required  or  dissipation  must  be  added  to  the  model.  We  have  tested  the 
two  approaches  and  demonstrated  that  adding  a  dissipative  term  to  the  model  leads  to  a 
significant  damping  of  polymer  stress  at  small  scales.  Fig.  5  Since  the  drag  reducing  effect 
is  significant  at  scales  not  much  larger  than  the  smallest  flow  scale,  an  upwind  high-order 
scheme  combined  with  local  artificial  dissipation  was  adopted  (Dubief  et  al,  2004a). 

DNS  in  channel  flows  allowed  the  statistical  characterization  of  LDR  (low  drag  reduction) 
and  HDR  (high  drag  reduction)  regimes,  as  shown  in  Figs.  6  and  7  (Dubief  et  ai.,  2004&). 

As  expected  from  experiments  (Warholic  et  aZ.,  1999),  DR=35%  produces  a  velocity 
profile  whose  log  region  is  only  shifted  upward,  while  DR=47%  and  DR=60%  show  significant 
changes  in  the  slope  of  the  log  law.  Although  we  did  not  attempt  to  determine  an  accurate 
limiting  value,  the  transition  between  LDR  and  HDR  seems  to  be  near  40%  of  drag  reduction, 
in  agreement  with  Warholic  et  al.’s  measurements,  The  turbulent  velocity  fluctuations  in 
the  transverse  directions  exhibit  a  different  behavior  than  in  the  streamwise  direction.  The 
peak  of  (a'  denotes  the  root  mean  square  of  a  fluctuations)  shifts  away  from  the  wall  and 
its  magnitude  increases  slowly  compared  to  the  Newtonian  flow  when  normalized  by  Ur»  The 
wall-normal  component  Uy  follows  an  opposite  trend;  behaves  as  Uy  and  consequently  does 
not  appear  on  the  plot  for  clarity.  In  c^ag  reduced  flow,  the  maximum  of  is  higher  or 
comparable  with  the  DR=0%  case,  as  found  in  experiments  (Warholic  et  aZ.,  1999;  Ptasinski 
et  alj  2003;  White  et  al.^  20036).  The  strong  reduction  of  the  transverse  fluctuations  suggests 
that  polymers  target  preferentially  the  vortices,  as  they  produce  significant  fluctuations  of 
Uy  and  Uz, 

The  balance  of  shear  stress, 

-  ^1  -  — ^  +  (1  -  P)Txy  =  0  ,  (5) 

contains  significant  information  regarding  the  mechanism  of  polymer  interaction  with  the 
mean  flow.  In  Eq.  (5),  Txy  is  the  average  of  As  indicated  in  Fig.  7,  the  Reynolds  shear 
stress  is  reduced  with  increasing  drag  reduction.  For  DR=60%,  is  approximately  a 

third  of  its  magnitude  in  the  Newtonian  case.  Conversely,  the  polymer  stress  increases  with 
increasing  drag  reduction.  At  DR=60%,  its  near-wall  contribution  to  Eq.  (5)  has  the  same 
magnitude  as  the  Reynolds  shear  stress,  thus  showing  that  polymers  have  a  significant  mean 
effect  in  the  mechanism  of  HDR  flows,  as  argued  by  Warholic  et  al  (1999). 

To  characterize  the  energy  which  is  stored  or  released  by  polymers  from  each  fluctuating 
velocity  component,  we  consider  the  Reynolds  stress  equation  (where  summation  does  not 
apply  to  the  subscript  a) 

^dtu\  =  {-UadjUaUj)  -{-{-Uada)p-^  -^UadjdjUa  +  ^-^UadjTaj  .  (6) 

Aa  Pc.  '  ^  ^ 

The  relation  between  polymers  and  turbulence  enters  through  the  polymer  work  Bay  product 
of  the  velocity  and  polymer  body  force,  which  could  dampen  <  0)  or  enhance  (E^  >  0) 
the  energy  carried  by  velocity  fluctuations  Uq..  The  other  terms  denote  work  by  advection 
Aa,  pressure  redistribution  Pa  and  viscous  stress  Va\  the  sum  of  these  terms  will  be  later 
referred  as  to  Newtonian  work.  > 

By  plotting  the  joint  probability  density  functions  of  polymer  work  with  other  works, 
we  were  able  to  correlated  negative  polymer  work  (polymer  storing  turbulent  energy)  with 
upwash  and  downwash  flows  and  positive  polymer  (polymer  giving  energy  back  to  the  flow) 
work  with  high  speed  streaks  very  close  to  the  wall  (see  Fig.  8).  The  latter  phenomenon 


6 


was  found  to  occur  in  regions  where  viscous  work  dominates  and  tends  to  damp  streamwise 
velocity  fluctuations.  By  releasing  their  energy  in  these  regions,  polymers  actually  overcome 
viscous  work  and  generates  a  growth  of  streamwise  velocity  fluctuations  (Fig.  8).  This 
transfer  of  energy  from  polymers  to  turbulence  is  responsible  for  preventing  the  flow  from 
relaminarizing  therefore  plays  a  key  role  in  MDR  (maximum  drag  reduction). 

3.9  DNS  of  boundary  layer 

The  simulation  of  turbulent  boundary  Jayer  started  late  in  phase  I  and  extends  to  phase  II. 
Statistical  results  for  the  velocity  (not  shown  here)  are  consistent  with  channel  flow  simu¬ 
lation.  Here  we  present  the  coupled  case  with  higher  viscoelasticity  which  was  thoroughly 
investigated  for  its  interesting  properties.  After  the  stationary  state  in  time  was  reached, 
averaging  was  performed  over  approximately  220  inertial  time-units  (4  flow-through  times). 
The  predicted  drag  reduction  along  the  streamwise  length  of  the  plate  exhibits  initially  a 
slight  drag  increase  close  to  the  inflow  boundary,  which  evolves  rapidly  to  a  larger  drag  re¬ 
duction,  reaching  40%  close  to  the  outflow  location,  where  a  small  region  of  constant  drag 
reduction  is  observed  (Fig.  9).  In  addition,  the  conformation  shows  that  there  is  a  larger 
region  close  to  the  inflow  where  stretching  is  most  intense  and  that  this  causes  more  rapid 
damping  of  the  vortices  downstream.  This  results  in  weaker  stretch  downstream,  which 
reaches  a  constant  value  close  to  the  outflow  of  the  current  computational  domain  (Fig.  10). 
In  addition,  the  concept  of  polymer  stretch  being  clearly  out  of  phase  with  its  result  on  the 
flow  field  (for  example,  vortex  intensity)  is  reinforced  (Fig.  11).  The  evolution  of  the  stretch 
of  the  polymer  in  these  coupled  simulations  underlines  the  importance  of  direct  numerical 
simulations  of  viscoelastic  turbulent  boundary  layers  because  they  capture  the  experimen¬ 
tally  observed  development  region  (White  et  a?.,  2003a)  even  for  the  case  where  the  polymer 
is  modeled  with  assumption  of  homogeneous  concentration.  Channel  flow  simulations  can 
only  capture  a  spatially  stationary  region,  which  does  not  provide  a  complete  understanding 
of  the  evolution  of  drag  and  microstructure  that  can  be  applied  toward  optimization  of  an 
injection  experiment. 

This  work,  its  validation  and  preliminary  results  are  available  in  an  internal  report  Dim- 
itropoulos  et  al  (2003). 

3.10  Drag  reduction  with  fibres 

Drag  reductions  of  up  to  21%  were  calculated,  with  the  largest  drag  reductions  observed  for 
non-Brownian  rods  at  semi-dilute  concentrations.  These  findings  suggest  that  elasticity  is 
not  necessary  in  order  to  achieve  Results  and  fibre  parameters  are  reported  in  Paschkewitz 
et  al  (2004). 

The  mean  velocity  and  RMS  velocity  fluctuations  for  the  Newtonian  and  drag-reduced 
cases  are  presented  in  figure  12.  The  fibre  drag-reduced  mean  velocity  profile  was  found 
to  be  similar  to  that  observed  in  other  drag-reduced  flows,  as  for  instance  polymer  flows: 
the  slope  of  the  log-law  region  is  roughly  the  same  as  that  of  the  Newtonian  case,  but  the 
intercept  with  the  linear  region  (u“^  =  y+)  is  shifted  upwards.  This  log-law  shift  is  indicative 
of  a  thickening  of  the  viscous  sublayer.  The  RMS  velocity  fluctuations  are  also  qualitatively 
similar  to  those  seen  in  polymeric  simulations,  with  an  increase  and  outward  shift  of  the 
streamwise  intensity  and  a  decrease  in  the  wall-normal  and  spanwise  fluctuations. 

To  determine  the  details  of  the  fibre-flow  interactions  suggested  by  the  scatter  plots  and 
single  point  correlations,  we  consider  instantaneous  visualizations  of  the  flow  field  and  the 
correlated  fibre  quantities.  In  figure  13  we  present  a  projection  of  the  velocity  vectors  and 
contours  of  normal  fibre  stress  in  the  plane  perpendicular  to  the  flow  on  the  left,  and  on 
the  right  we  show  contours  of  the  fibre  wall-normal  force  directly  affects  the  flow  at  the 
same  location  along  with  the  velocity.  Fibre  stress  is  large  near  a  stagnation  region  between 
vortices,  and  a  dipole  of  fibre  body  force  results  from  this  monopole  of  stress.  The  region  of 
positive  fibre  y-body  force  is  associated  with  a  region  in  which  the  velocity  vectors  point  in 
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the  negative  y-direction  and  vice-versa  and  thus  acts  to  decelerate  the  surrounding  vortices. 
We  observe  similar  behavior  for  T33  as  shown  in  figure  14.  Again,  the  largest  value  of  stress 
is  associated  with  a  stagnation  region  between  vortices  and  gives  rise  to  body  forces  that 
oppose  the  vortex  motion.  These  decelerations  in  turn  weaken  vortex  structure  and  disrupt 
the  near- wall  vortex  generation  mechanism  responsible  for  skin  friction. 

4  Discussion 

This  work  has  led  to  physical  explanations  of  the  mechanism  of  polymer  and  fibre  drag 
reduction  using  a  multiscale  approach  in  the  polymer  study. 

The  joint  pdfs  presented  in  Fig.  8a  and  8b  shows  how  polymers  enhance  the  streamwise 
momentum  in  high  speed  streaks  located  around  y"*"  =  5.  Polymers  are  allowed  to  coil  in 
these  regions  due  to  the  reduction  of  turbulent  kinetic  energy  by  viscous  dissipation  as  shown 
in  Fig.  8b.  Extraction  of  energy  from  near-wall  vortices  by  polymers  occur  as  polymers  are 
pulled  around  the  vortices,  either  by  upwash  or  downwash  flows.  Using  Brownian  dynamics 
simulations,  Terrapon  et  al  (2004)  showed  that  polymers  experience  significant  straining 
around  vortices,  leading  to  large  stretching.  The  importance  of  each  transfer  of  energy  for 
the  self-sustained  turbulence  in  drag  reduced  flows  has  been  established  by  simple  numerical 
experiments  reported  in  Dubief  et  al  (2004a), 

Based  on  these  results,  we  propose  to  include  the  effects  of  polymers  in  the  autonomous 
regeneration  cycle  put  forward  by  Jimenez  &  Pinelli  (1999).  This  cycle,  shown  in  Fig.  15, 
explains  how  wall  turbulence  is  self-sustained  through  mean  shear,  non-linear  interactions, 
near-wall  vortices  and  streaks,  in  the  Newtonian  case.  Polymers  fit  at  the  center  of  this  cycle 
by  extracting  energy  from  the  vortices  and  releasing  energy  in  the  streaks.  The  stretching 
of  polymers  is  governed  by  the  mean  shear  and  non-linear  interactions  as  shown  in  Terrapon 
et  al  (2004).  This  simple  mechanism  appears  to  apply  to  LDR  and  HDR  regime,  as  a 
matter  of  fact,  it  should  apply  to  any  regime  where  streaks  and  vortices  are  present.  Fig. 
15  also  contains  a  sketch  of  the  different  actions  of  polymers  around  a  vortex.  Only  the 
most  dramatic  phenomena,  vortex  damping  and  near-wall  high-speed  streaks  enhancement 
are  highlighted. 

For  fibres,  instantaneous  visualizations  of  the  flow  field  (not  shown)  show  that  high  stress 
regions  are  confined  to  intervortex  regions  and  are  most  pronounced  in  regions  of  intense 
vortex  activity.  These  regions  disappear  when  the  vortex  structures  are  weakened.  The  high 
stresses  were  found  to  be  caused  by  strong  planar  alignment  of  fibres.  We  therefore  infer 
that  the  fibres  are  reoriented  in  the  flow  direction  after  the  vortex  structures  are  weakened. 
Visualizations  of  the  principal  axis  of  the  fibre  orientation  tensor  as  a  function  of  time 
demonstrated  that  this  is  the  case.  Some  high  stress  regions  persist  even  after  the  vortices 
are  weakened  due  to  the  lag  in  this  reorientation  process.  After  the  fibres  reorient  in  the 
flow  direction,  the  reduction  in  local  fibre  stress  allows  the  vortices  to  reemerge  and  the 
turbulence  is  sustained  in  a  weakened  state. 

5  Important  Finds  and  Conclusion 

The  present  simulations  were  able  to  provide  extremely  valuable  insights  in  pol3aner  dynamics 
and  in  the  transfer  of  energy  between  polymers  and  turbulence.  The  numerical  method  was 
shown  to  apply  to  different  geometries  (boundary  layer)  and  other  additives  (fibres).  The 
most  important  results  which  led  to  the  models  of  polymer  and  fibre  drag  reduction  are 
listed  here: 

1.  At  low  We  rare  events  of  strong  biaxial  fiows  cause  the  chain  to  unravel  (Terrapon 
et  al.,  2004). 

2.  at  higher  We  different  flow  types  can  stretch  the  polymer  molecule,  but  only  bursts  of 
biaxial  extensional  flow  can  fuUy  extend  it  (Terrapon  et  al,  2004). 
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3.  Simulations  of  LDR  and  HDR  regimes  were  performed(Dubief  et  al^  20046). 

4.  Transfer  of  energy  between  polymers  and  turbulence  were  found  to  occur  in  a  well- 
organized  manner  triggered  by  coherent  turbulent  structures(Dubief  et  al,  20046). 

5.  Drag  reducing  and  turbulence  enhancing  eflFects  caused  by  polymer  dynamics  were 
shown  to  be  confined  to  the  near-wall  region  (Dubief  et  a/.,  20046). 

6.  The  mechanism  of  polymer  drag  reduction  has  been  described  as  a  modified  auto¬ 
regeneration  cycle  of  near  wall  turbulence  (Dubief  et  al,  20046). 

7.  A  new  algorithm  for  the  simulation  of  viscoelastic  turbulent  boundary  layer  flow  was 
developed  to  simulate  friction  drag  reduction  of  incoming  Newtonian  turbulence  for  a 
homogeneous  polymer  solution  (Dimitropoulos  et  al.y  2003). 

8.  major  viscoelastic  effects  are  located  with  20%  of  the  boundary  layer  thickness  (Dim¬ 
itropoulos  et  ai.,  2003).  ^ 

9.  A  phase  difference  between  polymer  stretch  and  vortex  damping  was  found  to  exist 
and  be  dependent  upon  the  polymer  relaxation  time  (Dimitropoulos  et  a/.,  2003). 

10.  Drag  reduction  with  fibre  additives  has  the  same  effect  on  statistics  as  polymer  addi¬ 
tives  at  LDR.(Paschkewitz  et  aZ.,  2004). 

11.  For  fibres,  elasticity  has  a  negative  impact  on  drag  reduction  effectiveness(Paschkewitz 
ef  flZ.,  2004). 

12.  Fibres  create  drag  reduction  via  an  extensional  mechanism  in  inter-vortex  regions  (Paschkewitz 
et  al,  2004). 

13.  The  drag  reduction  effectiveness  of  rigid  fibres  is  limited  due  to  requirement  of  the  in¬ 
teraction  of  multiple  vortices  for  drag  reduction  to  occur.  As  drag  reduction  decreases, 
vortices  are  weaker  and  move  farther  apart  (Paschkewitz  et  aZ.,  2004). 

6  Significant  Software  Development 

Three  codes  were  developed  and  validated  during  this  work:  a  channel  flow  code  coupled  with 
a  Brownian  dynamics  solver,  a  channel  flow  solver  for  constitutive  polymer  or  fibres  models 
and  a  boundary  layer  code  for  constitutive  polymer  models.  All  codes  are  fully-parallelized. 

7  Special  Conditions  for  Further  Research 

The  constitutive  polymer  model  used  (FENE-P)  did  not  incorporate  internal  modes.  Since 
polymers  are  not  significantly  stretched  at  HDR,  the  lack  of  internal  modes  which  will  in¬ 
crease  polymer  stress  compared  to  the  FENE-P  prediction  is  very  is  very  likely  to  be  respon¬ 
sible  for  an  underestimation  of  the  magnitude  of  drag  reduction.  Further  effort  should  be 
made  toward  the  development  of  a  more  complete  macroscopic  model. 

Fibres  hold  the  hope  of  combining  with  polymers  to  produce  higher  drag  reduction  as  it 
was  suggested  by  some  experiments.  The  coupling  of  the  two  additives’  effects  is  currently 
investigated. 

The  physical  insights  gained  by  our  numerical  simulations  can  be  used  to  derive  or  assess 
new  macroscopic  models  for  DNS,  Large  Eddy  Simulation  (LES)  or  Reynolds  Average  Navier 
Stokes  (RANS)  simulations. 
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Table  1:  Polymer  parameters  used  for  the  viscoelastic  simulations.  The  Weissenberg  number 
Wcro  is  normalized  the  wall-shear  stress  for  the  Newtonian  simulation  (DR=0%), 


Figure  1:  Square  extension  as  a  function  of  the  distance  to  the  wall  for  the  four  different 
models.  Closed  symbol:  mean,  open  symbol:  rms.  a  :  Constitutive  FENE-P,  □  :  Brownian 
FENE)-P,  <  ;  FENE,  v  :  Multichain  (5  modes)  . 


Figure  2:  Comparison  of  the  different  models  along  the  same  trajectory,  o  :  a  :  Consti¬ 

tutive  FENE-P,  o  :  Brownian  FENE}-P,  <  :  FENE,  v  :  Multichain. 
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(d)  (e)  (f) 


Figure  3:  Conditional  average  of  a*,  Q  and  R  for  the  polymer  molecules  crossing  the  thresh¬ 
old  value  q/qmax  =  r  as  a  function  of  the  time  At  before  crossing  this  threshold.  {a,b,c): 

r  =  0.65  and - ,We  =  1.0;  — ,  We  =  1.5; - ,  We  =  3.O.;  ■■■■,  We  =  6.0.  (d,e,/): 

We  =  3.0  and - ,r  —  0.95; - ,  r  =  0.85; - ,  r  =  0.75.;  — ,  r  =  0.65.  Shown  for 

comparison  as  a  thin  continuous  line:  (a-c)  results  from  a  bead-spring  chain  with  Nt  =  fl 
beads  at  We  =  1.5;  (d-f)  results  from  the  FENEl-P  model  for  r  =  0.95. 


13 


Figure  4:  Instantaneous  view  of  the  lower  half  of  the  channel  showing  the  isosurface  Q  =  1,9 
(grey)  representing  the  vortices,  the  isosurface  (turquoise)  of  cr*  =  1.6  and  the  polymer 
molecules  (red)  with  qjqmax  >  0.8  at  We  =  3. 


Figure  5:  Spectra  of  polymer  stress.  Treatment  of  the  advection  term  in  the  FENE-P: 
upwind  scheme  combined  with  local  artificial  dissipation  (blue);  addition  of  a  dissipative 
term  to  the  model  (red). 
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Figure  6:  Left:  Mean  velocity  profiles  scaled  with  inner  variables.  - :  MDR  asymptote, 

U'^  =  11. 7 In  +  17.  Symbols  and  lines  are  defined  in  Table  1.  Right:  RMS  of  velocity 

fluctuations  scaled  with  inner  variables..  Newtonian  simulation  (DR=0%):  -  , 

- ,  .  For  the  viscoelastic  simulations,  symbols  are  defined  in  Table  1,  is  denoted 

by  symbols  only;  is  indicated  by  symbols  connected  by - . 


Figure  7:  Reynolds  shear  stress  ( — --  ,  o  and  ■  as  in  Tab,  1)  and  polymer  stress  ( . 

DR=35%, - :  DR=60%)  normalized  by  Ur  and  i/,  DR=47%  is  not  plotted  for  clarity. 
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Figure  8;  Joint  probability  density  functions  of  polymer  work  vs.  other  relevant  quantities. 
For  each  quantity,  fluctuations  are  normalized  by  the  respective  standard  deviation.  Only 

data  at  DR=60%  is  shown.  Fig.  a;  Ex  vs.  Ux  ( - )  and  Uy  ( - )  at  j/+  =  5.  Fig. 

b:  Ex  vs.  Nx  =  Ax +  Px  +  K  ( - )  and  Idtul  ( - )  at  =  5.  Fig.  c:  Ey  vs.  Ux 

( - )  and  Uy  ( - )  at  y+  =  50.  Fig.  d:  Ey  vs.  Ny  =  Ay  +  Py  +V„  ( - -  )  and  ^dtul 

( - )  at  J/+  =  50. 


16 


Figure  9;  Coupled  DNS  for  p  =  0.9,  Ree  =  300,  =  3600  and  Wcr  =  50:  Predicted  drag 

reduction  along  the  streamwise  length  of  the  boundary  layer. 
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Figure  10:  Coupled  DNS  for  /?  =  0.9,  Ree  =  300,  =  3600  and  Wcr  =  50:  Average 

polymer  extension  along  the  wall  normal  direction  at  different  locations  downstream. 
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Figure  11:  Coupled  DNS  for  13  =  0.9,  Ree  =  300,  =  3600  and  Wtr  =  50:  Instanta¬ 

neous  snapshot  of  the  isosurfaces  of  the  second  invariant  of  the  velocity  gradient  tensor  with 
Q=0.075  (grey)  and  of  the  trace  of  the  conformation  tensor  with  threshold  1600  (black). 


Figure  12:  Mean  velocity  {left)  and  RMS  velocity  fluctuations  (right)  scaled  with  friction 
velocity.  Fibre  parameters:  Pe=1000,  aspect  ratio=100,  nL^  =  287 
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Figure  13:  Contours  of  fibre  stress  T22  (left)  and  fibre  y-body  force  (right)  with  y-z  plane 
projection  of  velocity  vectors 
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Figure  14:  Contours  of  fibre  stress  T33  (left)  and  fibre  z-body  force  (right)  with  y-z  plane 
projection  of  velocity  vectors 
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Figure  15:  Left:  Sketch  of  the  cycle  of  wall  turbulence  regeneration  with  energy  transfer  from 

the  polymers  to  the  flow  and  vice  versa, . denotes  the  action  of  polymers  on  turbulence, 

- the  main  actions  of  stretching.  Right:  Vortex  pumping  fluid  from  the  near-wall  region 

and  creating  turbulence-damping  polymer  work  and  re-injecting  stretched  polymers  into  the 
near-wall  region,  thereby  generating  turbulence  enhancing  polymer  work. 
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